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Exact  models  of  finite  end-item  population,  finite  repair  capaci¬ 
ty  repairable-item  systems  are  developed  using  Markov  process  analyses 
for  both  transient  and  steady  state  environments.  Unlike  most  currently 
used  multi-echelon  models,  the  infinite  population,  infinite  repair  ca¬ 
pacity  restrictions  are  removed.  Exponential  failure  and  repair  times 
are  assumed  and  the  system  is  modeled  as  a  closed  Markovian  queuing 
network. 

In  the  transient  case,  the  finite  set  of  differential  equations, 
and  in  the  steady-state  case,  the  finite  set  of  difference  equations, 
are  solved  by  numerical  techniques.  The  adequacy  of  these  techniques 
for  yielding  solutions  to  practical  systems  is  also  discussed. 
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1.  Introduction 

Consider  a  typical  multi-echelon  repairable-item  inventory  system 
as  shown  schematically  in  Figure  1.  Shown  there  is  a  two  location 
(bases) ,  two  level  of  supply  (spares  at  bases  and  depot) ,  two  level  of 
repair  (base  and  depot)  system  which  we  shall  denote  as  a  (2,2,2)  sys¬ 
tem.  The  nodes  BUi  (i  =  1,2)  represent  operating  and  spare  units 
(we  consider  for  now  only  a  single  item  such  as  a  final  assembly  or  a 
key  component)  at  base  i  ,  BRi  (i  =  1,2)  represent  the  repair  fa¬ 
cility  at  base  i  ,  DU  represents  depot  spares,  and  DR  the  depot 
repair  facility. 

Our  goal  is  to  develop  exact  mathematical  models  for  such  finite 
calling  population  (finite  number  of  items),  finite  repair  capacity, 
repairable  item  provisioning  systems  in  both  time-varying  and  steady- 
state  environments.  Specifically,  we  wish  to  find  the  state  probability 
vector  (the  probability  distribution  for  the  system  being  in  its  various 
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Figure  1.  Multi-echelon,  repair¬ 
able  item  system. 

possible  states)  which  will  allow  us  to  then  calculate  measures  of  per¬ 
formance  such  as  availability  (the  probability  that  at  least  some  desir¬ 
able,  prespecified  number  of  components  is  operational) .  Ultimately, 
these  models  will  be  used  to  yield  the  optimal  combination  of  spares  and 
repair  channels  at  each  location  in  the  system. 

Assuming  times  to  component  failure  and  component  repair  times 
to  be  exponentially  distributed  random  variables,  we  have  a  continuous 
time  Markov  process  (CTMP) .  The  process  is  driven  by  a  rate  matrix 
Q  =  {^ij}  5  where  is  the  "rate"  of  going  from  state  i  to  state 

j  ;  that  is,  letting  X(t)  represent  the  system  state  at  time  t 

Pr{X(t+At)  =  ,j  1 X ( t )  =  iTf  .  ,  . 

At  5  i  r  J  ; 

q .  .  =  -  £  q  .  .  . 

11  4  ^ij 

3  J 

(i*j) 

For  example,  suppose  the  (2,2,2)  system  pictured  in  Figure  1 
is  in  a  state  (call  it  i  )  for  which  the  depot  spares  pool  is  not 
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empty  (at  least  one  spare  is  on  hand  at  the  depot) .  Suppose  we  consider 
the  event:  a  component  fails  at  base  1.  Describing  this  state  i  by 

the  vector  (nBul>  nBR;[y  nBU2’  nBR2 ’  ^U’  ’  where  denotes  the 

number  of  components  at  node  k  in  the  "network,"  this  event  takes  the 
system  to  a  state  j  ,  namely,  (n^,  nBR1,  nBU2,  nBR2>  =,,,,-1,  n^+l)  , 
at  the  rate  =  Aa^  n^^  ,  where  1/A  is  the  mean  time  to  failure 

(MTTF)  of  a  component  and  is  the  probability  (or  percentage)  of 

failed  items  requiring  depot  repair. 

If  we  denote  the  state  probability  (row)  vector  at  time  t  by 
Tr(t)  -  (jr^(t),  TT^  ( t )  ,  .  ..,  TTg(t)3  9  that  is,  the  ith  element,  TT^  ( t)  , 
is  the  probability  of  the  system  being  in  state  i  at  time  t  (there 
is  a  finite  number  of  states  [call  this  number  S]  even  though  this  number 
can  be  quite  large),  then  we  must  solve  the  finite  set  of  first-order, 
linear  differential  (Kolmogorov)  equations 

tt*  (t)  =  TT(t)Q  .  (1) 

For  steady-state  solutions,  we  are  required  to  solve  the  finite  set  of 
linear  algebraic  steady  state  equations, 

o  =  TTQ  ,  (2) 

where  tt  =  (tt^,  .  ..,  tt^)  is  the  steady-state  probability  vector  and 

0  is  a  row  vector  of  all  zeroes.  In  both  steady-state  and  transient 
cases  we  have  the  further  condition  that  the  probabilities  sum  to  one, 
namely, 

1  =  £  (t)£  =  US  9 

where  e  is  a  column  vector,  with  all  components  equal  to  1. 
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2.  Transient  Environment 


We  are  often  interested  in  what  happens  to  such  systems  in  a  time- 
varying  environment.  For  example,  a  sudden  increase  in  effort  (say  a 
peacetime  to  wartime  shift)  may  cause  a  sudden  decrease  in  MTTF.  In 
such  situations,  it  is  necessary  to  have  7T ( t)  ,  and  we  must  solve  the 
finite  set  of  first  order  linear  differential  equations  given  in  (1) . 

Except  for  very  small  systems  (one  or  two  states)  analytical  tech¬ 
niques  such  as  Laplace  transforms  are  intractable.  Since  we  have  a  fi¬ 
nite  set  of  equations,  numerical  methods  can  be  employed.  Numerical  in¬ 
tegration  schemes  such  as  Runge-Kutta  or  predictor-corrector  methods 
are  possibilities.  We  choose  a  different  approach,  however,  which  is 
referred  to  by  some  as  randomization ,  and  has  been  shown  to  be  more  ef¬ 
ficient  for  these  kinds  of  problems  [see  Arsham,  Balana,  and  Gross  (1983) 
or  Grassmann  (1977a)].  For  details  on  this  technique,  which  can  be 
derived  by  a  probabilistic  argument  when  viewing  the  CTMP  in  a  certain  way, 
see  Grassmann  (1977a  and  b)  or  Gross  and  Miller  (1984a  and  b) . 

The  computational  formulas  are  as  follows.  Consider  a  discrete 
time  Markov  chain  (DTMC)  with  single-step  transition  probability  matrix 

P  =  Q/A  +  I  , 

where 


A  =  max  q . . 

1  xx1 

l 


that  is,  A  is  the  maximum  of  the  absolute  values  on  the  diagonal  of 
the  Q  matrix.  Since  a  diagonal  element  of  Q  is  the  negative  of  the 
sum  of  the  other  elements  in  the  row  (rows  of  the  Q  matrix  sum  to 
zero),  A  is  actually  the  absolute  value  of  the  minimum  (largest 
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negative)  diagonal  element  of  the  matrix.  This  DTMC  is  referred  to  as 

fk) 

a  uniformized  embedded  DTMC  of  the  CTMP .  Denoting  by  <J>  the 

state  probability  vector  of  this  DTMC  after  k  transitions,  it  can  be 
shown  (see  the  above  cited  references)  that 


TTj  (t) 


l 

=n  J 


k 

kT 


-At 

e 


For  computational  purposes,  it  is  necessary  to  truncate  the  infinite 
sum.  The  truncation  error  can  be  easily  bounded  since  we  are 
discarding  a  Poisson  "tail,”  so  that  the  computational  formula 
becomes 


where 


ITj  (t) 


T(t,e) 

l 


k=0 


(k)  (At)k  e“At 
J  k! 


(3) 


T(t ,e)  =  min  N: 


N 

1 

n=0 


e~At(At)n 

n ! 


>  1  -  £ 


e  being  the  maximum  tolerable  error  (specified  by  the  user) .  One 
advantage  of  this  method  over  numerical  integration  is  an  exact  bound 
on  the  computational  error. 

The  major  computational  effort  in  using  (3)  is  now  reduced  to 
finding  the  state  probability  vector,  ?  0f  the  uniformized  em¬ 

bedded  DTMC.  This  can  be  readily  accomplished  by  the  usual  recursion, 


4(0)  -  Tr(O)  ;  .  4<»p  . 


(4) 


Gross  and  Miller  (1984a)  give  a  more  efficient  procedure  than  the  suc¬ 
cessive  vector-matrix  multiplication  of  (4) ,  which  takes  advantage  of 
the  sparsity  of  the  P  matrix. 
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3.  Steady-state  Environment 

Solving  for  the  steady-state  probability  vector  TT  requires 
solving  the  set  of  linear  algebraic  equations  of  (2) .  Since  one  of  these 
equations  is  redundant,  it  is  necessary  to  reduce  the  equation  set  by 
one  and  use  1  =  TTe  as  the  final  equation.  Thus  (2)  can  be  reformulated  as 

b  =  TTA  ,  (5) 

where  b  is  a  vector  of  all  zeroes,  except  for  the  last  element,  which 
is  a  1,  and  A  is  the  Q  matrix  with  the  last  column  replaced  by  l's. 

For  relatively  small  systems,  the  solution  can  be  obtained  by 
inverting  A  to  get 

TT  =  bA_1  . 

However,  for  most  realistic  problems,  the  state  space  (and  hence  dimen- 
sion  of  the  A  matrix)  is  too  large  to  obtain  A-1  efficiently  or 
accurately.  This  situation  suggests  iterative  procedures  such  as 
Jacoby  or  Gauss-Seidel . 

Consider  the  A  matrix  as  a  sum, 

A  =  L  +  D  +  U  , 

where  L  is  a  lower  triangular  matrix,  D  is  a  matrix  with  only  diagonal  ele¬ 
ments,  and  U  is  an  upper  triangular  matrix.  Then  (5)  can  be  written  as 

tt(L+D+U)  =  b 
or 

ttD  =  b  -  tt(L+U)  .  (6) 

We  can  use  (6)  in  an  iterative  fashion, 

7r(n+1)D  =  b  -  T[(n)  (L+U)  , 

-  6  - 
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where  we  begin  the  procedure  with  some  initial  guess,  say  7T 


(0) 


This 


TT 


procedure  is  called  Jacoby  iteration.  Note  that  in  performing  the  cal¬ 
culations,  since  D  is  a  diagonal  matrix,  we  compute  ,  7T^n+^  , 

2 ,  ...  successively.  If,  as  we  compute  the  ,  we  replace 

the  7T^n^  on  the  right-hand  side  [e.g.,  in  computing  7T^n"*"'^  >  the 
Tf^n^  vector  is  modified  to  be  tt^  =  ,  ^l^1^  9  ■  5 

^j+l5  •••>  9  this  procedure  is  referred  to  as  Gauss-Seidel 

iteration,  and  in  matrix  representation  is 


(UT+D)  Jn+1')  =  b  -  LT7r^  , 

T  T 

where  tt  and  b  are  now  column  vectors,  and  U  ,  L  are  the  trans¬ 


(8) 


poses  of  U  and  L  ,  respectively. 

Two  questions  remain  to  be  answered  concerning  use  of  the  itera¬ 
tive  procedures  of  (7)  or  (8);  namely,  (i)  do  the  procedures  converge, 
and  (ii)  when  should  the  iterations  be  terminated?  In  general,  these 
procedures  may  not  necessarily  converge,  although  for  our  well-structured 
Markov  process  convergence  will  take  place.  The  stopping  criterion  gen¬ 
erally  used  is  the  Cauchy  criterion,  namely,  stop  when 


max 

i 


TT(n+1)  -  TTU) 
i  i 


<  £ 


0  5 


(9) 


where  £q  is  an  "arbitrarily"  chosen  small  number.  We  found  using  the 
fractional  difference  version  of  (9),  namely,  stop  when 


max 

i 


TrCn+l)  _  (n) 

i  i 


.(n) 


<  £, 


(10) 


to  be  somewhat  more  effective.  While  there  has  been  some  success  in 
using  Gauss-Seidel  (G-S)  on  Markov  models  [see  Kaufman,  Gopinath,  and 
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Wunderlich  (1981) ] ,  problems  exist  with  respect  to  rate  of  convergence 
and  appropriate  stopping  criteria.  The  G-S  convergence  rate  can  often  be 
improved  by  using  overrelaxation,  that  is,  by  weighting  with  a  coeffi¬ 
cient  greater  than  one  the  'n'in+1\  frfn+1\  ...,  used  in  calcu- 

0  1  j-1 

lating  TT^"^  [see  Kaufman,  et  at.  (1981)  or  Maron  (1982)]. 

Usually,  the  G-S  procedure  is  applied  to  a  set  of  equations  with 
a  nonsingular  matrix  (such  as  A  ) .  Consider  a  nonsingular  matrix  M 
with  positive  diagonal  elements  and  negative  off-diagonal  elements. 

The  G-S  procedure  is  known  to  converge  for  sets  of  equations  with  such 
an  M  matrix  [see  Varga  (1963)].  Now  consider  equation  set  (2),  namely, 

0  =  ttQ  . 

Multiplying  through  by  -1  gives 

0  =  tt  [  — Q  ]  , 

where  -Q  has  positive  diagonal  elements  and  negative  off-diagonal  ele¬ 
ments.  However,  it  is  singular,  since  one  equation  of  this  set  is  re¬ 
dundant.  Suppose  we  arbitrarily  set  tts  (assuming  there  are  S  states) 
to  one,  remove  the  last  row  of  the  Q  matrix  (call  this  reduced  matrix 
Q  ) ,  and  consider  solving  the  reduced  S-l  x  S-l  set  of  equations 

0  =  tt[-Q]  . 

A 

Now  -Q  is  an  M  matrix  and  convergence  is  guaranteed.  Of  course  the 
resulting  tk  values  are  relative  to  IT  =  1  so  that  they  must  be  renor- 

s 

malized  by  dividing  each  by  TT.  •  How  fast  convergence  takes  place 

still  is  a  key  question,  however.  It  turns  out  [see  Kaufman,  et  at. 
(1981)]  that  working  with  the  full  Q  matrix,  even  though  it  is  singu¬ 
lar,  speeds  convergence,  and  this  is  what  we  also  do. 


8 


T-490 


Another  procedure  is  to  use  the  uniformized  embedded  DTMC  of  the 
randomization  procedure  with  transition  probability  matrix  P  =  Q/A  +  I  . 
This  Markov  chain  has  limiting  probabilities  given  by 

4  =  &p  >  (ID 

and  they  are  identical  to  the  tt  of  the  CTMP  we  seek  [<J>  =  j>P  =>  = 

<j>C(Q/A)+l3  =>  0  =  £(Q/A)  =>  0  =  =  0  =  ttQ  ]  .  Solving  the  set  of 

equations  given  by  (11)  is  no  easier,  of  course,  than  solving  that  of 
(5) .  However,  we  know  from  Markov  chain  theory  that  limiting  probabil¬ 
ities  of  a  DTMC  can  be  found  by  iteration,  namely, 


TT 


(n+1) 


(12) 


Here  again,  we  have  computational  problems  associated  with  iteration, 
but  we  know  from  Markov  chain  theory  that  convergence  is  guaranteed  due  to 
the  existence  of  a  steady  state  vector  TT  (the  P  matrix  is  irreducible) . 
The  problem  of  when  to  stop  the  iterations  remains,  however.  Using  the 
Cauchy  criterion  here  results  in  problems  similar  to  those  found  when 
using  it  for  G-S  iteration,  namely,  successive  probabilities  can  differ 
by  very  small  amounts  and  still  be  far  from  the  steady  state  values. 

Wallace  and  Rosenberg  (1966)  provide  a  considerably  better  stop¬ 
ping  criterion  than  the  Cauchy  criterion  of  (10)  .  Their  stopping 
rule  is  based  on  estimating  the  rate  of  convergence  by  estimating 
the  second  eigenvalue  of  P  ,  and  turns  out  to  be:  "Stop  when 


(n+1)  (n) 

JL ~  TT  . 

f 

1  - 

\ 

(n+1) 

(n> 

> 

'  l/n 

TT 

-TT 

J 

For  details  of  this  development,  see  Wallace  and  Rosenberg  (1966)  or 
Gross,  Kioussis,  Miller,  and  Soland  (1984). 
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4.  Results 

The  following  section  gives  a  brief  summary  of  results  to  date. 
For  greater  detail,  we  refer  the  reader  to  Gross  and  Miller  (1984b)  and 
Gross,  Kioussis,  and  Miller  (1984)  for  the  transient  case  and  to  Gross, 
Kioussis,  Miller,  and  Soland  (1984)  for  the  steady-state  case. 

4.1  Transient  Case 

The  largest  system  solved  to  date  using  equation  (3)  directly 
was  a  (2,2,2)  system  (as  pictured  in  Figure  1)  with  18  components  at 
base  1  (of  which  4  were  spares)  ,  13  at  base  2  (of  which  3  were  spares)  , 
and  3  spares  at  the  depot.  The  base  repair  shops  had  2  parallel  service 
channels  each,  and  the  depot  repair  facility  had  four.  This  gave  a 
state  space  of  20,748  (Q  =  20,748  x  20,748). 

The  time-varying  environment  scenario  is  shown  in  Figure  2.  At 
time  6  ,  a  shift  in  MTTF  (1/A)  occurs  but  it  takes  until  time  10  for 
the  repair  facilities  to  Mcatch  up”  in  MTTR  (l/]i)  .  This  simulates  a 
change  in  usage  due  to ,  say,  a  shift  from  peacetime  to  wartime.  The 
measure  of  effectiveness  calculated  is  the  availability  at  time  t 
(t  =  1,2,..., 15)  ,  where  availability  is  defined  as  follows: 

A^(t)  E  Pr{at  least  14  components  are  operational  at  base  1  at 
time  t } * 

A^(t)  E  Pr{at  least  10  components  are  operational  at  base  2  at 
time  t } 

A.^2 ( t)  =  Pr{at  least  14  components  at  base  1  and  at  least  10 

components  at  base  2  are  simultaneously  operational 
at  time  t}  . 
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Figure  3  shows  a  plot  of  A^(t)  versus  t  .  Plots  of  A2(t)  and 
A^(t)  are  similar  in  nature.  The  graph  shows  an  initial  A^(0)  of  1.0 
(we  assume  at  time  zero  clIX  components  are  operational)  and  thereafter 
a  drop-off  toward  the  steady-state  availability  as  time  increases.  At 
time  6,  the  increase  in  failure  rate  occurs  and  A^(t)  begins  to  drop 
off  at  a  higher  rate,  heading  for  a  new,  lower  steady-state  availability. 
However,  the  increase  in  repair  rate  at  time  10  causes  A^(t)  to  begin 
to  rise,  heading  back  toward  the  original  steady-state  availability.  This 
run  took  approximately  25  minutes  of  CPU  time  on  a  VAX  11/780  computer 
using  the  randomization  computation  of  (3)  with  the  efficient  procedure 
given  in  Gross  and  Miller  (1984a)  for  calculating  t 

As  the  systems  become  more  complex  (more  bases,  multiple  component 
types,  indenture,  more  echelons,  etc.)  the  state-space  grows  rapidly.  We 
have  solved  a  problem  with  three  bases,  yielding  a  state-space  of  size 
43,278,703,  by  truncating  the  state-space  ("lumping"  low  probability 
states  into  several  absorbing  states  resulting  in  a  truncated  state  space 
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Figure  3.  A^(t)  versus  t  for  sample  run. 


of  approximately  15,000  states)  on  the  VAX  11/780  in  approximately  30 
minutes  [see  Gross,  Kioussis,  and  Miller  (1984)]. 
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4,2  Steady-state  Case 

Ironically,  computational  success  has  been  far  more  elusive  for 
the  steady-state  situation  than  for  the  transient  case.  The  problem  is 
the  stopping  criterion  for  these  iterative  procedures  (a  problem  not 
present  when  dealing  with  transient  solutions).  In  the  transient  case, 
the  randomization  procedure  guarantees  an  accuracy  to  within  a  pre¬ 
specified  £  .  For  steady  state,  using  either  the  Cauchy  or  the 
Wallace-Rosenberg  stopping  rule  does  not  guarantee  errors  within  £  . 
Table  1  shows  some  computations  for  a  (1,1,1)  system  which  is  the 
classical  machine  repair  with  spares  model  of  queueing  theory.  For  this 
model  the  availability  can  be  computed  analytically,  which  allowed 
us  to  estimate  the  actual  error.  The  columns  under  P-WR  show  the  results 
of  using  (12)  with  the  stopping  criterion  of  (13) ,  the  Wallace-Rosenberg 
approach,  while  the  GS-C  columns  show  results  for  (8)  with  -the  stopping 
criterion  of  (10) ,  the  Gauss-Seidel  approach. 

The  circled  elements  show  the  cases  for  which  the  error  specifi¬ 
cation,  £  ,  was  exceeded.  While  there  were  more  cases  of  exceeding  the 
stopping  rule  error  specification  in  P-WR,  the  error  excesses  were 
larger,  especially  for  the  larger  population  cases,  under  GS-C.  But 
GS-C  stopped  in  far  fewer  iterations  in  almost  all  cases  (except  for  the 
very  small  population  cases) ,  and  it  is  the  number  of  iterations  that 
consumes  most  of  the  CPU  time. 

The  last  column  shows  a  rerun  of  GS-C,  ignoring  the  stopping  cri¬ 
terion  and  performing  the  same  number  of  iterations  as  used  for  the  P-WR 
procedure.  The  errors  essentially  went  to  zero,  which  indicates  that 
if  a  better  stopping  criterion  could  be  found,  Gauss-Seidel  iteration 


13  - 


RESULTS  FOR  (1,1,1)  MODEL 


T-490 


co 

<r 

rH 

co 

CT\ 

n 

vO 

r^ 

o 

On 

o 

rH 

n 

rH 

n 

CN 

CO 

uo 

CO 

rH 

nt 

CO 

i — ( 

CT\ 

CO 

rH 

O 

CT\ 

CN 

rH 

ON 

in 

r^ 

co 

CN 

o 

vO 

o 

G> 

o 

CN 

nt 

i— i 

m 

On 

o 

CO 

00 

o 

co 

On 

O 

G\ 

vO 

P 
Ml  CJ 

Q.  r< 

S 


M  co 
•H  t — I 
ctf  0) 

a  e 
a)  a  3 

Cd  w 
=S=  CJ 


n 

o 

n 

n 

C3 

m 

m 

O 

n 

n 

O 

n 

n 

a 

n 

o  j 

rH 

rH 

o 

*H 

T«l 

d 

rH 

o 

-H 

d 

rH 

rH- 

Ci 

■— | 

rH 

co 

n 

co 

CO 

n 

CO 

a> 


CO 


m 


<r 


m 


o 

CN 


=fe  CO 

c 

T3  *H 
CU  4J  ^ 

•H  H  ^ 
co  a> 
a)  a 
P  o 


in 


o  o  o 

t — i  cn  in 


o 

o 

rH 


o 

o 


14 


^Percentage  of  time  the  population  is  operating  at  full  strength 
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might  be  a  viable  approach.  Runs  for  some  (2,2,2)  systems  and  more 
detailed  discussion  of  these  steady-state  procedures  can  be  found  in 
Gross,  Kioussis,  Miller,  and  Soland  (1984). 
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To  cope  with  the  expanding  technology,  our  society  must 
be  assured  of  a  continuing  supply  of  rigorously  trained 
and  educated  engineers.  The  School  of  Engineering  and 
Applied  Science  is  completely  committed  to  this  ob¬ 
jective. 


